/*
 heat-spectral lib. solves time-dependent PDEs using SEM-NI method.
 Copyright (C) 2009, Altieri Michele and Stefanoni Andrea
 
 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef EXACT_HPP
#define EXACT_HPP

#include<cmath>
#include<spectral>

// ******************************************************************
// Exact Problem
// ******************************************************************
// u_t - \Delta u = 0		in \Omega = [0,2\pi]^2 \times [0,T]
// u(x,y,0) = sin(2x) sin(2y)	in \Omega
// u(x,y,t) = 0			on \Gamma_D
//
// u(x,y,t) = exp(-8t) sin(2x) sin(2y)
// ******************************************************************

#define EXACT_PROBLEM ExactBilinear<spectral::Mesh>, Functional<spectral::Mesh,ExactFunc>


struct ExactSolution {
    double operator()(double x,double y,double t) const {
        return std::exp(-8*t)*std::sin(2*x)*std::sin(2*y);
    }
};

struct ExactSolutionX {
    double operator()(double x,double y,double t) const {
        return 2*std::exp(-8*t)*std::cos(2*x)*std::sin(2*y);
    }
};

struct ExactSolutionY {
    double operator()(double x,double y,double t) const {
        return 2*std::exp(-8*t)*std::sin(2*x)*std::cos(2*y);
    }
};

struct ExactFunc {
    double operator()(double x,double y) const {
	return 0;
    }
};

template<class MESH>
class ExactBilinear {
    public:
	ExactBilinear(MESH const & me): me(me) {
	    for(int i=0;i<4;++i) {
		bcond[i]=spectral::dirichlet;
	    }
	}

	template<class U,class V> inline
	    double operator()(U const & u,V const & v,int n) const {
		return
		    + me.prodL2(me.Dx(u,n),me.Dx(v,n),n)
		    + me.prodL2(me.Dy(u,n),me.Dy(v,n),n);
	    }

	// Initial data
	struct initial{
	    double operator()(double const & x, double const & y) const {
		return std::sin(2*x)*std::sin(2*y);
	    }
	};

	// Boundary conditions
	spectral::bctype bcond[4];

	// Function for the Dirichlet boundary condition
	struct dirichlet{
	    double operator()(double const & x, double const & y) const {
		return 0;
	    }
	};

	// Function for the Neumann boundary condition
	struct neumann{
	    double operator()(double const & x, double const & y) const {
		return 0;
	    }
	};

	// Methods on
	template<unsigned int N, typename T=void>
	struct on : public dirichlet {};

    private:
	MESH const & me;
};

const spectral::Point<double> LA(0.0,0.0);
const spectral::Point<double> LB(2.0*M_PI,2.0*M_PI);

#endif
